%% import data from .out files
% 
% label='1114';
% all_wind=importdata(['all_Lx37erg_z0.01_wind_',label,'.out']);
% all_rb=importdata(['all_Lx37erg_z0.01_rb_',label,'.out']);
% fir_wind=importdata(['first_Lx37erg_z0.01_wind_',label,'.out']);
% fir_rb=importdata(['first_Lx37erg_z0.01_rb_',label,'.out']);

% load('Lx37_z0.02_all.mat');
% 1  2   3  4   5  6   7   8   9  10 11 12  13    14 15 16  17  18 19
% m1,m2,Lx,ecc0,a,tb0,mt1,mt2,kw,kw2,mx,mx2,eccx,tbx,f,nnn,uuu,m_tr,i
% m1,m2,Lx,ecc0,a,tb0,mt1,mt2,kw,kw2,mx,mx2,eccx,tbx,f,m_transferrate,b_iso,m_tr,i
% 20 \delta T




%% initial parameters
global tmax ge AllmassNUM M1MAX M1MIN...
    CHAbeamingFlag SFR qNUM aNUM minLx maxLx alpha eddfac
alpha=2.7;
tmax=1d4*1d6; % in units of year
AllmassNUM=1d3;
M1MAX=150;
M1MIN=2;
SFR=GetSFR(50,alpha,0.535,0,M1MIN,5);
% SFR=10.6984;
% qNUM=10;
aNUM=100;
yearsc=3.1557d+07;
Msun=1.9891d30;
minLx=39;
maxLx=45;
% xi=Mdot/(sum(sum(data(:,1:2)))/tmax)/tmax*1d6;
CHAbeamingFlag=0;
eddfac=1d4;
ge=30;
%% delete some useless data

% [all_rb_1,all_wind_1]=DELallrbwind(all_rb,all_wind);



%% all in all. ONE for ALL!!!
allinall=[all_rb_1;all_wind_1];

global root filename
root=[pwd,'\result20201115_JOB20201114\'];
mkdir(root);


%% plot
for kw = 1 : 10
filename=['RocheLobekw_low',num2str(kw),'\'];

data=all_rb_1(all_rb_1(:,9)==kw,:);
plot_savedata(data);
end

for kw = 1 : 10
filename=['windkw_low',num2str(kw),'\'];

data=all_wind_1(all_wind_1(:,9)==kw,:);
plot_savedata(data);
end



for kw = 13:14
filename=['RocheLobekw_high',num2str(kw),'\'];

data=all_rb_1(all_rb_1(:,10)==kw,:);
plot_savedata(data);
end

for kw = 13:14
filename=['windkw_high',num2str(kw),'\'];

data=all_wind_1(all_wind_1(:,10)==kw,:);
plot_savedata(data);
end






function SETplot(picName)
global root filename
title(picName);
set(gca,'FontName','Times New Roman','FontSize',25,'LineWidth',1.5);
saveas(gcf,[root,filename,picName,'.bmp']);

end
function plot_savedata(data)
global root filename
if(isempty(data))
    fprintf('Empty')
    
    return;
end
mkdir([root,filename]);

plot_Lx_NNsum_gedian(data,1);
SETplot('Luminosity function integral')
ylim([0 2.2])

plot_Lx_N_gedian(data,1);
SETplot('Luminosity function differential')

plot_mass_tb_gedian(data,100,100);
SETplot('orbit and M2')
plot_beaming_gedian(data);
SETplot('beaming factor')
% save([root,filename,'data.mat'],'data');
end

